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ABSTRACT 



We consider the idealized expansion of an initially self-gravitating, static, 
singular, isothermal cloud core. For t > 0, the gas is ionized and heated to a 
higher uniform temperature by the formation of a luminous, but massless, star 
in its center. The approximation that the mass and gravity of the central star 
is negligible for the subsequent motion of the Hll region holds for distances r 
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much greater than ~ 100 AU and for the massive cloud cores that give rise to 
high-mass stars. If the initial ionization and heating is approximated to occur 
instantaneously at t = 0, then the subsequent flow (for r ^ 100 AU) caused by 
the resulting imbalance between self-gravity and thermal pressure is self-similar. 
Because of the steep density proflle (p oc r~^), pressure gradients produce a shock 
front that travels into the cloud, accelerating the gas to supersonic velocities 
in what has been called the "champagne phase." The expansion of the inner 
region at i > is connected to the outer envelope of the now ionized cloud core 
through this shock whose strength depends on the temperature of the Hll gas. In 
particular, we find a modified Larson-Penston (L-P) type of solution as part of the 
linear sequence of self-similar champagne outflows. The modiflcation involves the 
proper insertion of a shock and produces the right behavior at inflnity {v 0) 
for an outflow of flnitc duration, reconciling the long-standing conflict on the 
correct (inflow or outflow) interpretation for the original L-P solution. 

For realistic heating due to a massive young central star which ionizes and 
heats the gas to ~ 10^ K, we show that even the self- gravity of the ionized gas 
of the massive molecular cloud core can be neglected. We then study the self- 
similar solutions of the expansion of Hll regions embedded in molecular clouds 
characterized by more general power-law density distributions: p oc r~" with 
3/2 < n < 3. In these cases, the shock velocity is an increasing function of the 
exponent n, and diverges as n 3. We show that this happens because the 
model includes an origin, where the pressure driving the shock diverges because 
the enclosed heated mass is infinite. Our results imply that the continued photo- 
evaporation of massive reservoirs of neutral gas (e.g., surrounding disks and/or 
globules) nearby to the embedded ionizing source is required in order to maintain 
over a significant timescale the emission measure observed in champagne fiows. 

Subject headings: Hydrodynamics — Molecular Clouds — Hll Regions — Stars: 
Formation 



For a spherically symmetric molecular cloud core, initially at rest, the size rg of the 
region that can be ionized, is given by the standard formula (Stromgren 1939): 



1. 



Introduction 
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Eq. (|l]) assumes ionization equilibrium and the "on the spot" approximation. In the above, 
He is the electron density; is the ion density; 0:2 is the recombination coefficient to the 
second energy level of hydrogen; A^^, is the rate of ionizing photons from the star, assumed to 
be a constant; and rg is the radius below which all of the gas in the original cloud core ma; 
be considered to have fallen into the center (perhaps via a disk) to make a star of mass M^, 
If the virial velocity (thermal, turbulent, or magnetohydrodynamic) supporting the original 
(neutral) cloud core before star formation is denoted by ai, order of magnitude arguments 
yields ro ~ ri, the Bondi- Parker radius of this neutral gas. 

The square of the sound speed in the Hll gas is generally appreciably larger than af ; thus, 
the equivalent Bondi-Parker radius of the ionized gas, 

GM, 

will be considerably smaller than ri. 



r2 = (3) 



For typical numbers, M^, ~ 25 Mq, ai ~ 1 km s""*^, 02 — 10 km s""*^, we have ri ~ 10^ AU 
3> r2 ~ 10^ AU, with both ri and r2 much bigger than the physical radius of the star. Much 
interior to r2, the ionized gas will empty into the star (or more likely, into a disk if it has 
even a slight amount of angular momentum); whereas much exterior to r2, the gravitationally 
unbound Hll gas will expand outward, if it has not already reached pressure equilibrium with 
the surrounding cloud. Since rg 3> r2, we may henceforth ignore the gravitational field of the 
star on the flow of the Hll region beyond tq, although for purposes of making contact with 
earlier theoretical work, we shall begin by not ignoring the self-gravity of this gas. Since the 
material inside tq of the initial density profile should have fallen into the star, the observed 
presence of appreciable amounts of ionized gas at intermediate radii, ~ 10^ AU in typical 
ultracompact Hll regions, is awkward to explain. We defer until §6 the discussion of the 
special kinds of models that are probably required to explain ultracompact Hll regions. 

Assume now that the molecular cloud core initially had a power-law distribution of gas 
density that extends essentially to infinity: 

p(r) = iTr"". (4) 

If < 3/2, the ultraviolet radiation is trapped within a finite radius rs, and the Hll region 
is said to be "ionization bounded" (see Osterbrock 1989). U n > 3/2, the Hll region can 



^In the case when N^, oc t^, Newman & Axford (1968) found self-similar solutions for the expansion of an 
ionization bounded Hii region in a uniform HI cloud. 
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be either ionization bounded or "density bounded." In the latter finite output 

of ultraviolet radiation can ionize an infinite volume of gas beyond tq. The dividing line 
between being ionization bounded and density bounded arises when the density constant K 
equals a critical value Kcr'- 



(2n - 3)r2"-3 A^, 



.1/2 



47ra2 



(5) 



where Hi is the mean weight per particle of the ionized gas, ttlh is the hydrogen mass, and 
np = ne = p/2nimH- 

We wish to compare the value of K^r with the value K^, implied by the assumption that 
the power law initially extended inward from tq as well as outward, but that the part 
inward of ro has fallen into the center (perhaps via a disk) to make a star of mass M^,: 



(3 - n )M, 
47rro 



3— n 



(6) 



Taking the ratio of eq. (0) to eq. (Kl), we get 



(3 - n)M, 



a2 



{2n - 3)4nr^N^ 



1/2 



(7) 



For M, ~ 25 Mq, tq ^ 10^ AU, iV, 
23 (3-r2)/(2n-3)V2. 



10 



49 



02 



2.6 X 10-13 cj^s g-1^ KJK^ 



It is remarkable that factors of such disparate orders of magnitude as the dimensionless 
quantities in the two square brackets of eq. combine to give a ratio within two orders 
of unity. Nevertheless, since K^, represents a rough estimate of K and > Kcr, this 
calculation formally indicates that the Hll regions of 25 Mq (and lower mass) stars are likely 
to be ionization bounded, at least initially before any expansion occurs. However, if we 
assume that A^* scales roughly as (as indicated by the results of Vacca et al. 1996), the 
expression on the right-hand side scales as M"^, indicating that the Hll regions of the most 
massive O stars may be density bounded from the start, especially if such stars are born in 
regions with density gradients close to n = 3. They will then develop champagne flows as 
follows. 

When K K^: < K^r, the ionization front (IF) created by the idealized instantaneous 
appearance of a star at t = rapidly moves to infinity and establishes an isothermal structure 
with T ~ 10'* K. After the passage of the IF, the cloud remains out of mechanical balance 
and the pressure gradients will produce an expansion of the whole cloud. Due to the density 
gradient the inner regions expand faster than the outer regions and a shock travels through 
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the cloud, accelerating the gas to supersonic velocities. This is known as the "champagne 
phase" (e.g. Bodenheimer, Tenorio-Tagle & Yorke 1979). Franco, Tenorio-Tagle & Boden- 
heimer (1990; hereafter FTB) studied the evolution of Hll regions embedded in molecular 
clouds with steep density gradients. High spatial resolution infrared and radio recombination 
line observations toward several sources have found ionized gas accelerating away from the 
central source in the manner expected of champagne flow models (e.g. Garay et al. 1994; 
Keto et al. 1995; Lumsden & Hoare 1996; Lebron et al. 2001). Note that in several of the 
observed compact Hll regions (e.g., 29.96-0.02, G32.80+0.19B, G61.48+0.09B1) the inferred 
rate of ionizing photons imply excitation by central stars with masses > SOMq. 

Density profiles in massive molecular cores have also been extensively studied observa- 
tionally (e.g. Garay & Rodriguez 1990; CaseUi & Myers 1995; Van der Tak et al. 2000; 
Hatchell et al. 2000; for a review see Garay & Lizano 1999). Even though the environment 
is possibly clumpy on scales of tenths of pc, density profiles are well approximated by power 
laws with 1 ^ n ;S 2. Theoretical models of the formation of massive stars within dense and 
massive cores have assumed power law exponents in this range (Osorio et al. 1999; McKee 
& Tan 2002). Recently, Franco et al. (2000) have argued that radio continuum spectra of 
ultracompact Hll regions indicate initial density gradients with 2 n ^ 3. Clearly, more ob- 
servations with high spatial resolution are necessary to reliably establish the density profiles 
of the sites of massive star formation. 

The purpose of this paper is to study by similarity techniques the "champagne phase" of 
expansion of Hll regions with power law density distributions. In §2, we formulate the outfiow 
problem in the case of the singular isothermal sphere (SIS) that has p oc r^^, including the 
effect of self-gravity. Tsai & Hsu (1995) found the outflow analogue of the inside-out collapse 
solution (Shu 1977), but in which the SIS is sent into expansion by an outward propagating 
shock. In §3 we show that the Tsai & Hsu (1995) solution is actually the limit of a family 
of outflow solutions when (01/02)^ — > 1 from below. Furthermore, the outflow solution with 
the particular ratio of (01/02)^ = 0.75 corresponds to a piece of the time-reversed Larson- 
Pcnston (L-P) collapse solution (Larson 1969; Penston 1969), but with a shock inserted to 
obtain the correct asymptotic behavior for large distances (or early times). For realistic 
heating after the passage of an IF, i.e., for realistic values of (01/02)^ ^ 1, we show that 
the self-gravity of the Hll gas can be neglected. In §4 we extend our study to the evolution 
of champagne flows in the case of density distributions with power law exponents in the 
range 3/2 < n < 3, neglecting self-gravity. In §5 we flnd that the self-similar models have 
a shock propagating at constant velocity into the ionized gas, in good detailed agreement 
with the models of FTB. In particular, the shock velocity diverges as n — >^ 3. We show that 
this happens because the formal treatment extends the inner radius of the calculation to the 
origin. In such a treatment, the mass of driving Hll gas diverges when n > 3. We perform 
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more realistic calculations in such cases that cut holes in the gas distribution for r < Tq. 
In §6 we summarize our conclusions, and we discuss the implications of our results for the 
problem of ultracompact Hll regions. Finally, in the appendices we show that, for the scales 
relevant to molecular cloud cores, the isothermal assumptions for the gas and the shock are 
valid. 



2. Governing Equations 

Consider a star forming cloud core with the density profile of the SIS 

where ai is the sound speed of the cloud ai t < 0. Imagine that Bit t — 0, the central star 
turns on and heats the entire cloud core to a uniform temperature appreciably higher than 
it had originally (perhaps by the rapid propagation of an ionization front). Let a = 02 > oi 
be the sound speed corresponding to this new temperatTirc. In order to keep the initial gas 
density distribution unchanged from i = 0~ to t = 0"*", it is convenient to write 

where e = (ai/a)^ < 1. The original temperature was appropriate for a self-gravitating SIS, 
but the new higher temperature makes the same density distribution over-pressured. The 
gas expands and a shock front propagates through the cloud, setting it into a champagne 
outflow. 

The expansion of self-gravitating champagne flows is governed by the continuity equa- 
tion, 

where p is the gas density and u is the gas velocity, and by the momentum equation, 

du du _ dp GM 
dt dr p dr ' 

where G is the gravitational constant and M is the mass of the gas inside the radius r. An 
energy equation is not required since we we have assumed an isothermal equation of state 
P — a^p, with shock jumps also occurring isothermally (see Appendix B) . 

We further simplify the treatment by approximating the central star as a source of light 
and heat (to keep the gas ionized and warm), but not of mass, and we extend the range of 
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the flow calculation right into the origin. We justify this cavalier treatment of the complex 
situation interior to r = Tq by noting that the first thing all champagne solutions try to 
establish is a central region of quasi-uniform pressure (see §§3 and 4). All the available 
ionized mass interior to r = ro is homogenized by the high pressure in less than a sound 
crossing time (because supersonic motions are involved), ro/a. Thus, all detail about the 
structure of the region interior to r = tq is lost in a time < 5 x 10^ yr, which is quite short 
compared to the duration of interest for champagne flows. Apart from initial transients, 
then, the evolution of quasi-uniform Hll regions of high over-pressure expanding into exterior 
regions with power-law density profiles should follow the description given by the self-similar 
treatment of this paper. An exception to this statement arises if the region interior to 
r = ro is a continuous source of additional ionized gas (derived, say, from the continuous 
rather than instantaneous photoevaporation of a circumstellar disk and/or remnant neutral 
globules). We ignore this exception for the present, but we return to the possibility in §6 
when we comment on the problem of ultracompact Hll regions. 



2.1. Similarity assumption 



Following Shu (1977) we introduce the similarity variable 

r 



X 



at ' 



and we define the reduced density 



the reduced mass 



and the reduced velocity 



/ X a(x) 



M{r,t) = ^m(x), 

Gr 



u{r, t) = av{x). 



(12) 

(13) 

(14) 
(15) 



Substituting these expressions in eq. (jT^) and ([TTp, one obtains two coupled first order 
differential equations (ODEs) for the reduced density a and velocity v 



[{v-xr-i] 



1 da 
a dx 

dv 



a [X — v) 

X 



dx 



{x — v)a 



X 



(x-v), 
(x-v), 



(16) 
(17) 
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while the reduced mass is given by 

m = x'^a{x — v). (18) 

We identify t = as the initial instant and require t, x, m, a and v all to be positive. 
Then, the critical lines where the LHSs and RHSs of eq. ( [T6| ) and (|1^) vanish lie on the locus 



2 , , 

Vc = Xc - 1, ac = —. (19) 

Xc 

A special case with an analytical solution that passes smoothly through the critical line at 
Xc = 3 was cited by Whitworth & Summers (1985; eq. [3.5]): 

a = I v = lx. (20) 

In dimensional form, this is the model for the Einstein-De Sitter universe (in the non- 
relativistic limit where <^ so that we can ignore the contribution of the gas pressure in 
determining spacetime curvature), with a Hubble "constant" and cosmic mass density given 
respectively hj H = 2/(3t) and p = 3H^/{87cG) = l/{6nGt'^). 

In the Hll region problem, the reduced central density a(0) = ao will not usually be 
tuned to the specific value ao = 2/3 required to make a smooth crossing of the critical line. 
Even more, the pressure (and density) homogenization of the central regions cannot extend 
instantaneously to all space because of the existing density gradient at t = 0. Nevertheless, 
if we ignore the central star (or more accurately, if we put it back to fill the hole inside 
To), the behavior of the Einstein-De Sitter solution is generic for the inner regions of the 
Hll region in that the dimensional central density will decline with time as because of 
the expansion of the flow toward the exterior regions of lower pressure and density. Thus, 
for arbitrary values of ao, we must solve eq. (|T6|) and (0) subject to the inner BCs: 

a = ao and v = 0, at x = 0. (21) 

A series expansion near the origin (a regular singular point of eqns. ||I6| and |T^) now yields 



ao /2 \ 2 /r.f-.\ 
a = ao + — --aox + ..., (22) 



6 V3 

V = -X + — (- - ao] x^ + . . . . (23) 

3 45 V3 / ^ ' 

The heating of the cloud at t = introduces an imbalance between self-gravity and 
pressure that will induce the propagation of a shock and an outward subsonic flow of the 
entire system for t > 0. We assume isothermal jump conditions 

{Vu-Xs){Vd-Xs) = 1, — = {vu-Xsf, (24) 
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where the subscript u (d) indicates the value of the reduced velocity and density upstream 
(downstream) of the shock and Xg is the position (and velocity) of the shock in similarity 
coordinates. From eq. ( [T^ ) the position and velocity of the shock in physical space are 
Ts = XgCit and Us = drg/dt = axg- 

The outer solution to which this shock attaches is given by one of the "plus solutions" 
of Shu (1977), and it has the asymptotic behavior 

2e 2(1 -e) , , 

a — *• — , V , as a; — *• cxo, (25) 

x^ X 

with e < 1. 



2.2. Numerical procedure 

We start the numerical integration of the downstream flow at small x, using the second- 
order expansions (123) for a given value of ao- At each integration point, we apply the 
isothermal jump conditions and a provisional upstream solution is obtained. A consistent 
solution is found when the upstream solution, integrated to large x (~ 10'^), fulfills the 
asymptotic BCs (^). Thus, for a given ao there is a corresponding unique value of e such 
that the full solution satisfies the BCs (|2T|) and (pSl) and the jump conditions, eq. (pi). 



3. Outflow Solutions for the SIS 

By varying the value of a^, we obtain a family of outflow solutions. For increasing 
values of ao, the shock velocity and shock strength decrease, and the shock attaches to outer 
solutions characterized by increasing values of e < 1. Figure |l| shows the reduced velocity v 
and density a versus x for ao = 2/3, 1, 1.67, 4, and 7.90, corresponding to e = 0.545, 0.640, 
0.755, 0.916, and 1. 



3.1. The Tsai & Hsu solution 

There is a maximum value ao = 7.90 which corresponds to the solution that attaches 
upstream to the static unperturbed SIS solution with e = 1 (thick solid lines in Figure [1|). 
In this case the shock moves at Ug = 1.34a, setting the unperturbed SIS into motion. Since 
e = 1, the temperature of the cloud has not changed at t = 0. Thus, the shock is driven by 
an infinitesimal pressure pulse at the origin and acquires finite strength as it races down a 
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density gradient. This self-similar solution, found first by Tsai & Hsu (1995), is the outflow 
counterpart of the "expansion wave solution" (Shu 1977) for the collapse of the SIS, but now 
the SIS is sent into expansion by an outwardly propagating shock. Thus, the SIS is unstable 
not only to self-similar "inside-out collapse" but also to "inside-out expansion" (Tsai & Hsu 
1995). (Of course, physically, there is no reason to expect a perturbational heating of the 
cloud core without having, first, inward collapse to form a central star.) This special outflow 
solution is the limit of such solutions as e — 1. 



3.2. The Larson-Penston solution 

For the particular value ao = 1.67 the flow downstream from the shockfront is coincident 
with a time reversed piece of the L-P solution originally proposed to describe infall. Without 
introducing a shock, the L-P solution crosses smoothly the critical lines at = 2.34 (dotted 
lines in Figure 0), the reduced density being tangent to the critical line, and achieves the 
asymptotic behavior 

8 85 

a — -- and v 3.28, for a; ^ oo. (26) 

X'^ 

The L-P solution is represented by dashed lines in Figure |l|. In contrast, our outflow solution 
has a shock at Xs = 1.80 and has the asymptotic behavior 

1.50 

a — > — — and f ^ 0, for x — > oo, (27) 

x^ 

which is less dense than the SIS at large x, as appropriate for an outflow solution. In fact, 
as discussed by Shu (1977) the time reversed L-P solution does not represent a proper wind 
solution, precisely because it is overdense at large radii with respect to the hydrostatic SIS. 



3.3. The Einstein-De Sitter solution 

For the special case ao = 2/3, the solution interior to the shockwave at Xg = 2.06 is the 
Newtonian analogue to an Einstein-De Sitter [Q = 1) universe. The pressure (and density, for 
an isothermal gas) homogenization of the interior is perfect in this case, and Hubble's law of 
expansion is also exactly satisfied before the gas flow reaches the radius r = Xgat. Upstream 
from the shock front, the flow corresponds to a breeze (see next subsection) blowing slowly 
through an e = 0.545 singular isothermal sphere. Notice that this use of the Einstein-De 
Sitter interior solution (connection to an exterior solution via a shockwave) differs from that 
shown in Fig. 1 of Hennebelle (2001) (connection via a weak discontinuity). 
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Unlike its cosmological application, the Einstein-De Sitter (inner) solution represents 
nothing special in the problem of the self-similar isothermal expansion of self-gravitating 
Hll regions. It merely demarcates the case (see Table 1) when the post-shock density 
switches from values lower than the central level ao (for ao > 2/3) to values higher than 
the central level (for < 2/3). This switch occurs with increasing heating of the Hll gas 
relative to its neutral precursor (smaller values of e = [01/02]^). With weak heating, running 
down a density gradient established by the self-gravitating equilibrium of the pre-ionization 
state and modified by a slowly blowing, upstream breeze counts for more than compression 
of this gas by a relatively weak shockwave. With strong heating, shockwave compression 
gains in the competition relative to farther propagation into regions of lower density. 



3.4. General Behavior 

With the exception of the "inside-out expansion" solution of Tsai & Hsu (1995), the 
upstream flows have Vu > 0. This is due to the mechanical imbalance following the heating 
of the cloud at t = causing the gas to expand subsonically before the shock arrives. Such 
subsonic behavior is characteristic of a "breeze" solution in the nomenclature of the solar- 
wind literature, and we have referred to it as such. We note, however, that the arrival of the 
shockwave pushes the breeze into supersonic expansion relative to the origin. The strength 
of the shock and the value of Vu depend on the new temperature of the cloud at t = 0"*". 
Table 1 summarizes the results for the self-similar models with different values of the central 
density Uq. The columns are: 1) the value of ao; 2) the parameter e of the upstream solution; 
3) the position of the shock, x^; 4) the post-shock velocity Vd, 5) the post-shock density a^; 
5) the pre-shock velocity Vu', and 6) the pre-shock density «„. The entries in this table are 
approximated to three significant figures. The rounding-off error affects the fulfillment of 
the jump conditions at the level of ^ 2%, a sufficient accuracy for practical applications. 



3.5. The limit e 

As discussed in § 0, when the cloud core is heated instantaneously at t = 0, the value 
of e must be chosen for t = 0^ so that the physical density p(r, 0) remains unchanged. For 
a cloud heated by the passage of an ionization front when a massive star turns on in the 
center of the core, one expects the equivalent thermal speed to increase from ai ~ 1 km s~^ 
to a ~ 10 km s~^. Thus, the value of e for the upstream flow should be quite small, e ^ 10^^. 
Since the equilibrium configuration of the cloud has e = 1 for t < 0, the condition e <^ 1 
implies that self-gravity is not very important in the subsequent evolution of the cloud for 
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t > 0. 

We find that as e — > 0, the reduced velocity v converges to an invariant solution with the 
fastest and strongest shock. Table 1 shows that most flow properties have converged pretty 
much to their limiting values when e ~ 10^^. The only formal exception is the reduced 
density, which converges to a shape invariant solution that scales as ao 0.47 e. Thus, in 
the limit e ^ 0, we find it convenient to re-scale the density as 

= (28) 
The physical density for t > in term of R is given now by 

Comparison with eq. shows that i? = 1 at t = O"*". 

The upper panel of Figure |^ shows the reduced velocity v in the limit e = 0. The 
position of the shock front is at Xg = 2.56 and the upstream reduced velocity, due to the 
cloud general expansion, is f„ = 1.02. The dotted line shows the locus of the critical line eq. 
(p!9|). The lower panel of Figure |^ shows the reduced density R. The dashed line shows the 
function R = S^x/xsY, which corresponds (see eq. [Q) to a region of uniform but steadily 
decreasing density, 

'^'') 

given by spreading the original gaseous mass interior to = Xgat evenly over the enclosed 
spherical volume. Except for a slight increase of R from its average interior value of 3 to the 
postshock value Rd = 3.20 just downstream from the shockfront at x = Xs, Figure |^ shows 
that the high gaseous pressure does a fairly good job of ironing out pressure differences in 
the interior volume. 

In summary, for a density distribution in mechanical equilibrium at t < 0, which comes 
out of balance at t = 0"*" because it has been heated suddenly to a temperature much 
higher than the original value, one expects that the effects of self-gravity can be neglected 
with respect to thermal pressure. In this limit where e — > 0, we verify explicitly that the 
self-similar solution converges to a unique solution, equivalent to the gravitationless case, 
G = 0. Thus, in the next section we will formulate the problem of self-similar champagne 
flow solutions for power laws of the gas density distribution different from n = 2 also in the 
limit of negligible self-gravity . 
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4. Generalized Formulation Ignoring Self-Gravity of the Ionized Gas 

In this section we discuss the champagne flows associated with generalized power-law 
density-proflles with p cx r~" in the limit where the self-gravity of the cloud can be neglected. 
For t < 0, these density proflles result from equilibrium between the self- gravity of the cloud 
core and other forces in addition to thermal pressure (see e.g. discussion on polytropic 
clouds in Galli et al. 1999). We assume that, after heating at t = 0, thermal pressure will 
dominate over such forces in the evolution of the cloud. The governing equations are given 
by eq. (11^) and (|lTD setting (7 = 0. We choose the same self-similar variable and reduced 
velocity deflned by eq. (|12[) and ([1^) but we generalize the density proflle of eq. as 

p{r,t) = ^R{x). (31) 



The nondimensional equations are 

di 



dR R 

[{v - x)2 - 1]— = -[v{v -x){n-2)-n], (32) 



[(v-xf-l]— = n. 33 

dx X 

For X ^ oo the boundary conditions are f — ^ and R ^ 1. Then, the asymptotic 
expansions are 

V - and R^l + — — ^-^ for x ^ oo. (34) 

At the origin, 

Tl 

V —X and i? oc x" for x — 0. (35) 
3 

This boundary condition on R implies that as r ^ the density is uniform and is only a 
function of time p oc This is the expected behavior of the central zone of the Hll region, 
where the sound crossing time is smaller than the expansion time. 



We apply the isothermal jump conditions (p^) to connect the upstream and downstream 
solutions. Since the equation for the reduced velocity is now decoupled from the equation for 
the reduced density, a simple way to flnd the position of the shock front, x^, is to integrate 
numerically eq. (^) outward from x = 0, and integrate inward from a large x. We check 
every point of the downstream solution until the jump condition in the velocity (P^ ) is 
fulflUed. 

Once the reduced velocity v is known, the density equation (|32|) can be integrated as 

J x' \ V - x') - 1 



Xl, 
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For the outer solution, = oo and -R^ = 1. For the inner solution, xi, = Xg and Rb = Rd, 
where Rd is the downstream reduced density evaluated from the jump condition (p^. 

Before proceeding to discuss the self-similar champagne flows, let us examine the case 
n = 3. In this case, eq. (|32D and (0) have the analytic solution R = Cx^ and v = x, where 
C is an arbitrary constant. The jump conditions eq. (|2^) imply that Xg ^ (yo as Vd —>■ x. 
Thus, in spatial coordinates the shock front, = x^at, and the shock velocity, Ug = aXs, go 
to infinity. 

In dimensional variables, the analytic solution for the case n = 3 simply corresponds to 
a Hubble flow in a = universe: u = r/t, p oc t~^. The propagation of the shockwave 
to infinity leaves behind it an expanding Hll region of uniform but ever decreasing density 
that looks like a Little Bang. We will discuss below the physical meaning (or lack of it) of 
this remarkable result. 



5. Outflow Solutions for Different Gradients 

The nongravitating case for n = 2 was shown in Figure |^. Figure |^ show the analogous 
reduced velocity v and density R for n = 2.99. The dotted line in the upper panel corresponds 
to the critical line. In this case the shock front is at Xs = 19.25. This number is also the 
velocity of the shockwave relative to the isothermal speed of sound of the Hll region. In 
contrast, the pre-shock velocity of the gas (relative to the origin) is only Vu = 0.16, because 
at a given spatial position there is less time for the cloud to expand before the shock arrives. 
The reduced density R in the lower panel has a correspondingly large post-shock density 
increase. The dashed line shows the function 

given by spreading the original gaseous mass interior to = Xgat evenly over the enclosed 
spherical volume. Figure |^ shows that the shock dynamics raises the immediate downstream 
value from the average expectation 3/(3 — n) = 300 ai x = Xs to the actual postshock value 
Rd = 368. 

Table 2 summarizes the results for the self-similar models with different exponents n in 
the density profile. The columns are: 1) the exponent n; 2) the position of the shock, Xs', 3) 
the post-shock velocity Vd] 4) the post-shock density Rd] 5) the pre-shock velocity f„; and 
6) the pre-shock density Ru- One can see that Xg increases as n — 3. 

Figure ^ shows the comparison of our self-similar solutions with the numerical results 
obtained by FTB in the case n = 1.7. We converted the values of the velocity and density 
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calculated by FTB at t = 1.72 x 10^ yr (see their Figure 3c) into our self-similar variables 
(eq. |T2|, |15| and |31|) using their values of the sound speed and their normalization for the 
density profile. Figure ^ shows the excellent agreement between their numerical simulation 
and our similarity calculation. 



5.1. Analytic Approximation 

Following FTB we compare the results of the self-similar champagne flows with a simple 
analytic approximation. Assume that the pressure driving the shock is uniform, 

P(t)=aV(t), (38) 
where p(t) is the mean density at a given time, i.e. 

Pit) = '° ' = p[rM with n < 3, (39) 







n 



where ^^(t) is the instantaneous position of the shock front. This must be equal to the 
(strong shock) post-shock pressure given by 

P{t) = ius-Uufp[r,it)], (40) 

where Uu is the upstream velocity at the position of the shock front. Equating these two 
pressures, one obtains 

x. = - = ^P^+-, (41) 
a \ 6 — n a 

which is equivalent to eq. (24) of FTB when Uu = a. Eq. ( PT| ) diverges as n — > 3 because, 
in this limit, the mass inside any radius diverges if the origin is included (eq. p9|) . 

Figure ^ shows \og{xs) vs. n for the self-similar models. The long dashed line corre- 
sponds to eq. ( ^TD where we have used m„ = a (as FTB). The dotted line instead is eq. (^Tj) 
with Uu = 0. The difference is due, in part, to the pressure driving the shock being not 
exactly uniform as was assumed in the approximation above. The deviation from uniform 
density was commented upon in §§3 and 5 when we discussed the approximate nature of the 
fits provided by the dashed lines in the lower panels of Figures ^ and |^. 

To avoid the divergence at the origin, one could integrate from the radius of influence tq 
(or include a uniform core as considered by FTB). In that case, the density can be evaluated 
for all n and is given by 

m = /[r.(t)]p[r,(t)], (42) 
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where 

/KW, 4 - M'"") (43) 

" \ iln|r,(()/r„f if n = 3 ^ ' 

and Q = 1 — [ro/rs(t)]^ is a volume factor. The shock velocity is then 

Us = a^f[rs{t)] + Uu. (44) 

For n < 3, /(r^) — > 3/(3 — n) when ^ tq; thus, the shock reaches the asymptotic constant 
value given by eq. (0). For the case n > 3, the shock accelerates to oo as ^ oo. The 
latter divergence arises only because a shockwave of finite outward momentum is allowed to 
run into spherical shells of ever decreasing mass. 

Finally, as discussed in Appendix A and B, even though the expansion velocities become 
very large as n 3, adiabatic cooling never offsets photoionization heating and the gas 
remains isothermal. We also show that for the scales of interest in molecular clouds, even if 
the shock velocity becomes large, the assumption of an isothermal shock remains valid. 



6. Discussion 

Compact Hll regions with velocity gradients indicative of champagne flows are character- 
ized by values of the emission measure (EM) in the range 10^ cm~^ pc ^ EM ^ 10^ cm~^ pc 
(e.g. Garay et al. 1994). The results of FTB and the self-similar models presented here show 
that when a molecular core is ionized and heated out of equilibrium, steep density gradients 
characteristic of star-forming regions produce shocks that travel at constant velocity and 
accelerate the ionized gas to supersonic speeds. This phase of supersonic expansion of the 
ionized gas poses a lifetime problem for champagne flows, more severe than the one pointed 
out by Wood & Churchwell (1989) in the case of ultracompact Hll regions. 

The emission measure of the self-similar models presented here at time t = 0^, just after 
ionization of the cloud core is 

EM„./;„.„A=(5;^)(5^)'rJ- (45) 

As we discussed in § |l|, if if < K^r the Hll region is density bounded and therefore able 
to develop a champagne flow. With K^r defined by eq. (^, and our fiducial values of the 
physical parameters for a 25 Mq star, this condition results in the following upper limit on 
the EM: 
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in the middle of the range indicated by the observations. However, inserting the expression 
of p{r,t) given by eq. ( pl| ) into the eq. (|i5|), one easily obtains EM(t) oc t^~^" for t ^ ro/a. 
Thus, the expansion of the champagne flows results in a rapid decrease of the emission 
measure, and in a short time the source fades away from observational classification as a 
compact Hll region. Inclusion of the effects of a fast stellar wind would only exacerbate the 
situation. 

These estimates show that a continuous source of ionized mass is required to keep 
up the density of the expanding champagne flow. Photoevaporation of circumstellar disks 
(Hollenbach et al. 1994; Richling & Yorke 1997) or mass loading by the photoevaporation 
and/or hydrodynamic ablation of remnant neutral globules surrounding the central star (e.g. 
Lizano et al. 1996; Redman et al. 1996) are natural solutions to maintain the high observed 
emission measures in champagne flows and, in general, in ultracompact Hll regions. In 
particular, proplyds in Orion, sites of low mass star formation, are known to have appreciable 
mass loss rates that can mass load the stellar wind (Storzer & Hollenbach 2000; Garcfa- 
Arredondo et al. 2002). Revealed O and B stars show a lower frequency of surrounding disks 
than stars of spectral type A and later (Natta et al. 2000; Lada et al. 2000), suggestive of 
a picture where the disks are photoevaporated during an earlier phase of evolution, perhaps 
as ultracompact Hll regions. As a first attempt to model quasi-spherical ultracompact 
Hll regions on this basis, it would be interesting to extend the methods of this paper to 
include a continuous source of ionized gas (due to photoevaporation of a disk or neutral 
globules into a mass-loaded stellar wind) that emanates at a prescribed rate and speed from 
the origin. 



7. Conclusions 

We have obtained self-similar champagne flow solutions for the expansion of power 
law gas density distributions after the gas has been uniformly heated out of mechanical 
equilibrium by the birth of a star at the center of a molecular cloud core. These solutions 
attach via a shock to upstream breeze solutions. 

In the case of the isothermal sphere with p oc r^^, the "inside-out expansion" found 
by Tsai & Hsu (1995) is the limit of the family of self-similar outflow solutions when the 
sound speed 02 of the Hll region is the same as the sound speed oi of the original molecular 
cloud core. The case (01/02)^ = 1 must include the effect of the self-gravity of the gas, 
and the outflow solution then attaches to the unperturbed static SIS upstream of the shock. 
Another member of this family, for (01/02)^ = 0.75, is a time reversed piece of the L-P 
collapse solution, with a shock allowing the upstream solution to have a correct outflow 
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( breeze) asymptotic behavior. 

For the high values of the gas temperature expected after the passage of an ionization 
front, (01/02)^ <IC 1, and the self-gravity of the gas can be neglected. The solution then 
approaches a shape invariant form. In the approximation that the self-gravity of the ionized 
gas can be ignored, we computed the self-similar champagne flows of Hll regions formed in 
molecular clouds characterized by power law density distributions with exponents 3/2 < n < 
3. These self-similar solutions behave as in the numerical models of FTB: in the "champagne 
phase" a shock moves with a constant speed into the ionized medium with a shock speed 
that increases with increasing density gradient. The speed of the shock diverges as n — 3 
because the mass of the driving gaseous piston diverges, if the origin is included. Instead, 
if the origin is excluded, the shock front velocity reaches an asymptotic constant value for 
n < 3. For n > 3 the shock accelerates to infinite velocity as Ts ^ 00, but only because finite 
outward momentum is inputted into spherical shells of ever decreasing mass. These results 
may help explain astrophysical champagne flows where expansion velocities are seen that 
are considerably larger than the sound speed 02 — 10 km s~^ associated with conventional 
Hll regions. (Driving by fast stellar winds may contribute to the perceived motions.) Despite 
the large expansion velocities produced in the case of the steepest pressure gradients (n 3), 
we show in the appendices that the isothermal assumption for the gas and for the shock are 
valid for the scales of interest in molecular clouds. 

The supersonic expansion of the ionized gas creates a severe lifetime problem for cham- 
pagne flows. A natural solution is the photoevaporation of circumstellar disks and/or rem- 
nant neutral globules which would help maintain the high observed emission measures in 
these sources. 

We thank Luis Rodriguez, Malcolm Walmsley and Frank Wilkin for useful comments on 
the manuscript. F. S. acknowledges support in the United States from the NSF and NASA, 
and in Taiwan from the National Science Council. S. L. and J. C. acknowledge support from 
DGAPA/UNAM and CONACyT. D. G. acknowledges support from grant COFIN-2000 and 
warm hospitality from the Instituto de Astronomia, UNAM in Morelia. 

A. Adiabatic Cooling vs. Photoionization Heating 

For an expanding nebula, the equilibrium temperature is given by the balance between 
photoionization heating Fph and radiative cooling Arad and adiabatic cooling Ag^. 
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We can estimate when adiabatic cooling becomes important by setting 

Tph = Aad. (Al) 
Assuming ionization equilibrium and the "on the spot" approximation, Fph is given by 



Tph ^ n^npa2-kT^, 



(A2) 



where is the stellar temperature (e.g. Osterbrock 1989). The rate of adiabatic cooling 
(for a constant mass) is given by 



A„A = — — = a 



V dt 



dp 



dt 



(A3) 



Then, the condition for adiabatic cooling to win over photoionization heating can be written 

as 



77-42^' 



(A4) 



where tree = l/'n.eQ!2 is the recombination timescale, and tp = p/\dp/dt\ is the timescale for 
the decrease in the density. 



Taking T^/T ~ 4, the critical ratio of timescales is 

^rec D 
tn 



(A5) 



On the other hand, using eq. ( ^I]) and (|39|) this ratio can be written as 



tp \ 3 J Ka2 



n-l 



(A6) 



Using the critical value for the ratio of timescales ( |A5D , one can solve this equation for 
a critical radius r^r beyond which the isothermal assumption breaks down. For convenience, 
one can write the density constant in eq. (|) as K = pmnnorQ, where p is the mean 
molecular weight and uq is the number density at the distance tq. Then, the critical radius 
can be written as 

ZjLJi^UjJbg \0 lb/ \lb/ 

Finally, for 0:2 — 2.6 x 10^^'^ cm'^ s~^ and the fiducial parameters ~ 10^ cm"'^ tq — 10^ 
AU, a ~ lOkms"^, and yU ~ 2, typical of massive molecular cloud cores, the critical radius 
becomes 

3 y/2 /3^ni/("-i) 
3 — 72/ \n, 



To 



780 



(A^ 
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where we used eq. ( ^T]) with = to obtain the last expression. 

This critical radius is large compared to tq and the sizes of molecular cloud cores where 
the stars are formed. Thus, adiabatic cooling never dominates over photoionization heating 
and the assumption of an isothermal gas is justified. 



B. Isothermal Shock 

We now examine the assumption that the shock is isothermal in the case when the shock 
speed is large (n ^ 3). The post shock temperature is 

^ , ^ X 10^ ( X K. (Bl) 



16 k ' VlOkms-i 

For large shock velocities the gas behind the shock will not be able to cool efficiently, thus 
the cooling time may become larger than the expansion time, tcooi > 4xp = rg/ug, and the 
shock may become energy conserving. 

The expansion time is 



,11 / ' \ I Ug ^ ^ 



The cooling time is given by 

_ 3 kTd 

tcool - 7 , [d6) 

where is the gas post-shock temperature, is the gas post-shock number density and 
Arad is the cooling function (in erg cm^ s^^). Substituting the post-shock temperature (|ET|), 
using eq. (|24D for the post-shock density and assuming that the upstream flow is as rest, the 
cooling time can be rewritten as 

9 i^f rriH a 

^cool 



32 /i rioArad Vro 

~ 1.4xl0^f -J^I^fL^^] f^i^Vs, (B4) 

V 10-23 erg cm'^ s^V \W A\J J ' ^ ' 

where we have again used the expression for K deflned in Appendix A. 

Thus, the ratio is 



^cool 7 ( "^^rcid ^ s \ / \ 

I 10-23 erg cm3 s'^ J UO^Auj UOkms-^ ' ' ^^^^ 
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This equation is expressed in terms of typical values of the cooling function (e.g. Dal- 
garno & McCray 1972) and shows that the ratio of cooling time to expansion time is small. 
Thus, the assumption of an isothermal shock will remain valid for the scales relevant in 
molecular clouds. 

REFERENCES 

Bodenheimer, P., Tenorio-Tagle, G., & Yorke, P. 1979, ApJ, 233, 85 

CaseUi, P., & Myers, P. 1995, ApJ, 446, 665 

Dalgarno, A., & McCray, R. A. 1972, ARA&A, 10, 375 

Franco, J., Kurtz, S., Hofner, P., & Testi, L. 2000, ApJ, L143 

Franco, J., Tenorio-Tagle, C, & Bodenheimer, P. 1990, ApJ, 349, 126 (FTB) 

GaUi, D., Lizano, S., Li., Z.-Y., Adams, F. C, & Shu, F. H. 1999, ApJ, 521, 630 

Garay, G., & Lizano, S. 1999, PASP, 111, 1049 

Garay, G., Lizano, S., & Gomez, Y. 1994, ApJ, 429, 268 

Garay G., & Rodriguez, L. F. 1990, ApJ, 362, 191 

Garcia- Arredondo, F., Arthur, J., & Henney, W. 2002, Rev. Mex. Astron. Astrof., 38, 51 

Hatchell, J., Fuller, G. A., Millar, T. J., Thompson. M. A., & Macdonald, G. H. 2000, A&A, 
357, 637 

Hennebelle, P. 2001, A&A, 378, 214 
Hunter, C. 1977, ApJ, 218, 834 

HoUenbach, D., Johnstone, D., Lizano, S., & Shu, F. H. 1994, ApJ, 428, 654 

Keto, E. R., Welch, W. J., Reid, M. J., & Ho, P. T. P. 1995, ApJ, 444, 765 

Lada, C. J., Muench, A. A., Haisch, K. E., Lada, E., Alves, J. F., ToUestrup, E. V., & 
Willner, S. P. 2000, AJ, 120, 3162 

Larson, R. B. 1969, MNRAS, 145, 271 

Lebron, M., Rodriguez, L. F., & Lizano, S. 2001, ApJ, 560, 806 



-22- 



Lizano, S., Canto, J., Garay, G., & HoUenbach, D. 1996, ApJ, 739 
Lumsden, S. L., & Hoare, M. G. 1996, ApJ, 464, 272 
McKee, C., & Tan, J. C. 2002, Nat. 416, 59 

Natta, A., Grinin, V. P., & Mannings, V. 2000, in Protostars and Plantes IV, eds. V. 
Mannings, A. P. Boss, & Russell, S. S. (Tucson: The university of Arizona Press), p. 
559 

Newman, R. C., & Axford, W. I. 1968, ApJ, 153, 595 
Osorio, M., Lizano, S., & D'Alessio, P. 1999, ApJ, 525 808 

Osterbrok, D. E. 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (Mill 
Valley: University Science Books) 

Penston, M. V. 1969, MNRAS, 144, 425 

Redman, M. P., Williams, R. J. R., & Dyson, J. E. 1996, MNRAS, 280, 661 

Richling, S., & Yorke, H.W. 1997, A&A, 327, 317 

Shu, F. 1977, ApJ, 214, 488 

Storzer, H. & HoUenbach, D. 1999, ApJ, 539, 751 

Stromgren, B. 1939, ApJ, 89, 526 

Tsai, J. C., & Hsu, J. J. L. 1995, ApJ, 448, 787 

Vacca, W. D., Garmany, C. D., & ShuU, J. M. 1996, ApJ, 460, 914 

van der Tak, F. F. S., van Dishoeck, E. F., Evans, N. J. II, & Blake, G. A. 2000, ApJ 537, 
283 

Whitworth, A., & Summers, D. 1985, MNRAS, 214, 1 

Wood, D. O. S., & Churchwell, E. 1989, ApJSS, 69, 831 

Yorke, H. W., Tenorio-Tagle, G., & Bodenheimer, P. 1983, A&A, 127, 313 



This preprint was prepared with the AAS IM^jX macros v4.0. 



-23- 



Table 1. SIS Outflow Solutions as Function of ao 
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Table 2. Flow Parameters as a Function of n 



n 












1.5 


2.15 


1.33 


1.85 


0.93 


1.25 


1.6 


2.23 


1.44 


2.04 


0.96 


1.28 


1.7 


2.30 


1.55 


2.26 


0.98 


1.30 


1.8 


2.38 


1.66 


2.52 


1.00 


1.32 


1.9 


2.47 


1.78 


2.83 


1.01 


1.33 


2.0 


2.56 


1.91 


3.20 


1.02 


1.34 


2.1 


2.66 


2.05 


3.65 


1.01 


1.35 


2.2 


2.77 


2.21 


4.21 


1.00 


1.34 


2.3 


2.91 


2.39 


4.93 


0.98 


1.33 


2.4 


3.07 


2.59 


5.87 


0.95 


1.31 


2.5 


3.27 


2.85 


7.17 


0.91 


1.29 


2.6 


3.55 


3.17 


9.10 


0.85 


1.25 


2.7 


3.96 


3.64 


12.26 


0.76 


1.20 


2.8 


4.67 


4.42 


18.48 


0.65 


1.14 


2.9 


6.33 


6.16 


36.94 


0.47 


1.08 


2.95 


8.76 


8.65 


73.73 


0.34 


1.04 


2.99 


19.25 


19.20 


367.57 


0.16 


1.01 



-25- 




1 2 3 4 5 

X 




1 2 3 4 5 



X 

Fig. 1. — Top panel: reduced velocity v for the outflow solutions of the SIS, including self- 
gravity. The different curves from top to bottom correspond to ao = 2/3, 1 , 1.67, 4 and 
7.90 {thick solid curve). The dotted line is the critical line v = x — 1, Bottom panel: reduced 
density a for the same values of a^. The dotted curve is the critical line a = 2/x. In both 
panels, the long- dashed curves show the L-P solution without a shock. The parameters of 
the solutions are described in Table 1. 
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Fig. 2. — Top panel: reduced velocity v for the exponent n = 2, in the hmit e = 0. The 
dotted line is the critical line v = x — 1. The shock front is at Xg = 2.56. The pre-shock 
velocity is Vi = 1.02. Bottom panel: reduced density R. The dashed curve corresponds to 
R = 3{x/xs)^ and shows the deviation from uniform density. 
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Fig. 3. — Top "panel: reduced velocity v for the exponent n = 2.99. The dotted line is 
the critical line v = x — 1. The shock front is at Xs = 19.25. The pre-shock velocity 
is only Vi = 0.16. Bottom panel: reduced density R. The dashed line corresponds to 
R = 300 (x/a;^)^'^^ and shows the deviation from uniform density. 
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Fig. 4. — Top panel: reduced velocity v for the exponent n = 1.7. The dotted line is the 
critical line v = x — 1. The shock front is at Xg = 2.30. The pre-shock velocity is = 0.98. 
Bottom panel: reduced density R. In both panels, the triangles show the results of the 
numerical models of FTB for n — 1.7 sit t — 1.72 x 10^ yr, expressed in our self-similar 
variables. 
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Fig. 5. — Plot of Xg vs. n showing the divergence of the shock position as n — > cxo. The dots 
represent the exact self-similar solutions. The long-dashed curve and the short-dashed curve 
correspond to the analytic approximation eq. ( ^T]) with Uu = a and Uu = 0, respectively. 



